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An effective analytical method for calculating energy spectra of cosmic-ray muons 
at large depths of homogeneous media is developed. The method allows to include 
an arbitrary (decreasing) muon spectrum at the medium boundary and the energy 
dependence of both discrete (radiative and photonuclear) and continuous (ionization) 
muon energy losses, with resonable requirements for the high-energy behavior of the 
initial spectrum and differential cross sections of the muon-matter interactions. 
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I. INTRODUCTION 



Cosmic-ray (CR) muons originate from the decay of unstable hadrons produced by the 
interactions of cosmic-ray primaries and secondaries with nuclei of the earth's atmosphere. 
Therefore the flux of CR muons contains information on primary cosmic rays (energy spec- 
trum, composition, anisotropy) as well as on some properties of particle interactions at high 
and super-high energies. 

During the last years the experimental investigations of CR muons with large low- 
background detectors for penetrating particles have expanded rapidly in a number of un- 
derground laboratories, in addition to the direct measurements in the atmosphere and at 
ground level. Side by side with the traditional range of problems of cosmic ray physics 
some additional aspects arise within the framework of investigations with the new facilities. 
Thus, for example, the flux of CR muons is used for calibration of the detectors and, at 
the same time, it is an important source of background events for the majority of under- 
ground experiments, especially in neutrino astronomy and astrophysics. Detailed study of 
this background is very important for further progress in astroparticle physics. 

Projects for the deep-underwater Cerenkov and acoustic detection of high-energy muons 
and neutrinos have been discussed for a long time. The ultimate aim of these projects is to 
build detectors of volume 10 7 — 10 9 m 3 or even larger which could be used, in particular, 
to detect muons of energy up to 10 3 TeV, in order to study the CR muon flux at energies 2 
to 3 orders of magnitude higher than those accessible in the present experiments. It should 
be particularly emphasized that detectors of so enormous volume can be used to accurately 
determine the energy, E, of individual muons passing through the apparatus if E exceeds a 
few TeV §. 

In the near future the initial stage of the underwater muon and neutrino telescopes 
DUMAND II in the ocean off the Hawaii island || and NT-200 in the Siberian lake Baikal 



j|] will be commissioned. Besides, several new programs have been proposed, such as a 
deep sea neutrino detector NESTOR in the southwest corner of the Peloponnisos || and 
the project AMANDA for a large scale muon and neutrino detector in deep ice at the South 
Pole || (see also Ref. @] for a short summary on next generation detectors). Precision 
calculations of different characteristics of the CR muon flux after propagation through thick 
water layers are an imperative element for the successful realization of these projects. 

The transport of high-energy CR muons through dense media has been the subject of 
theoretical investigations over many years with the use of analytical |8|-[16|, numerical []17|,[P8|1 , 



and Monte Carlo [|19| -|21| methods (see also Refs. p2|j23[| for a review of the early literature). 
In the majority of the papers listed the depth-intensity relation (DIR) was studied. However, 
for future experiments with large-scale underwater neutrino telescopes a detailed knowledge 
of the energy spectra of muons at very large depths (large zenith angles) will be required in 
addition to the total (oblique) intensities. Some results of calculations of the muon energy 
spectra at large depths of matter were presented in Refs. PJ20|,|2T] and in our previous 



papers |1J,|I5|, but the increasing requirements on accuracy of the calculations stimulate us 
to continue the investigation of the problem. 

The main difficulty in the calculation of muon intensity and spectrum at large depths 
consists in the fact that an ultrarelativistic muon of energy E above ~ 100 GeV can, with 
comparable probabilities, lose in a single event either a very small energy AE <C E or an 
energy AE ~ E with generation of a large electromagnetic or hadronic shower due to ra- 
diative or photonuclear interaction with matter. These fluctuations of the energy loss lead 
to a pronounced range straggling. Considering that the rate of radiative and photonuclear 
losses increases with energy, the fluctuation effect grows with energy and depth. An impor- 
tant consequence of this effect is the impossibility to define a threshold energy for a muon 
reaching a given depth. This fact presents a severe problem when reconstructing the surface 



(sea-level) muon spectrum from a measured DIR (see, e.g., Ref. 

Available exact analytical methods of cascade theory [II] require that the initial muon 



spectrum at the boundary of the medium be a power-law and that the differential cross 
sections for muon-matter interactions depend only on the fraction of energy lost by the muon, 
v = AE/E, but not on the muon energy, E, itself ("scaling"). It is also assumed that the rate 
of the continuous (ionization) energy loss is constant. We shall call this set of assumptions 
the SPS model (Scaling + Power-law Spectrum). Such model have been considered after 
Rozental' and Strel'tsov [J2£| in Refs. |8|-[T0|. Zatsepin and Mikhalchi || have suggested a 



very simple approximate solution to the muon transport equation (TE). Their approach has 
been generalized to a quasi-power initial spectrum [0]. The exact solution for DIR within 



the framework of the SPS model has been obtained by Nishimura |T0[ (see also Refs. |p"T|-|T3|j) 
with the use of the technique of integral transformations. Both these approaches have been 
employed with some modifications in numerous works (see, e.g., Refs. p~6|j2T[| ) . 

However, as is generally known, the assumptions of the SPS model hold roughly only at 
very high muon energies (above 1—10 TeV), and so the calculations based on the SPS solution 
should be corrected by one or another way. Despite the relatively weak energy dependence 
of the differential cross sections as well as the closeness of the real sea-level muon spectrum 
to a power-law form within wide energy intervals, these corrections prove to be very large 
and they increase with depth. The point is that the muon energy spectrum under thick 



1 In this connection the problem of prompt muons which appear in the atmosphere due to decays of 
charmed hadrons should be mentioned p5| , Pq1 . The data of the current underground experiments, 
those from European detectors (NUSEX, Frejus, MACRO), on the one hand, and those from 
Baksan and KGF, on the other hand, contradict each other (see Ref. f2(|). A certain part of these 
disagreements can be attributed to an inaccuracy in the computation of the fluctuation effect (see 
also Ref. 24] and a discussion in Ref. 



layers of matter depends exponentially on integrals of the differential cross sections with a 
weight which depends on the initial spectrum. To our knowledge there is not any proper 
and consistent way to calculate these corrections at large depths for the time being. 

In the present study we discuss a comparatively simple and universal method for the 
calculation of differential energy spectra as well as other important characteristics of CR 
muons at arbitrary depth, which allows us to avoid from the start the simplifying assump- 
tions about the scale invariance of the cross sections and the (quasi) power-law incident 
muon spectrum. The solution to the TE is constructed by iterations, starting from an 
initial approximation with the correct high-energy asymptotic behavior. In the range of 
applicability of the initial approximation (sufficiently high energies) it becomes feasible to 
introduce an (effective) analog of the threshold energy at the boundary which is very useful 
in many respects. One of the advantages of the computer realization of our approach (in 
comparison with the direct Monte Carlo simulation or a purely numerical technique) is its 
high performance which allows to carry out verifications of various hypotheses on the pri- 
mary CR spectrum and composition, charm production models, models of the photonuclear 
interaction, etc with good precision and in a negligible CPU time. This enables to estimate 
the sea-level muon spectrum using the data of underground/underwater measurements by 
exhaustion, avoiding to solve the much more difficult inverse scattering problem. 

It should be noted that the method under consideration is a development of our previous 



studies 14,15 



The organization of this paper is as follows. In Sec. [IT] we give some preliminaries and 
notations. We present also a very short review on some features of the differential cross 
sections for muon-matter interactions at high energies which will be needed later on. In 



Sec. [Ill] the solution to the TE in a continuous loss approximation is discussed briefly for 
the methodological goals. In Sec. fV| we consider the exact solution to the TE within the 
SPS model; for the present purpose (to study the asymptotic behavior of the TE solution 
at high energies) the simplest expression in the form of a series in powers of 1/E will be 
quite enough. In Sec. [V| we derive an approximate solution to the TE in the general case; 
essential properties of the solution are discussed in some detail and illustrated by the SPS 
model. The iteration algorithm for calculating corrections to the approximate solution is 
described in Sec. |VT] and the convergency of the algorithm is examined. Finally, in Sec. |V11 
we summarize the results and some perspectives for applications of the method. 



II. PRELIMINARIES 



A. Transport equation 

The propagation of relativistic muons through a homogeneous medium is described by 
the one-dimensional transport equation (TE) 

—D(E, x)- — [Pi(E)D(E, x)} = (D(E, x)) (2.1) 

with the boundary condition 

D(E,0)=D (E) . (2.2) 

Here D(E, x) is the differential energy spectrum of muons at depth x in the medium. In the 
general case 



x = sec# / p{z')dz' , 
Jo 



where p(z) is the density of the medium at distance z from the boundary, and $ is the 
angle of incidence measured from the normal to the boundary (zenith angle). The function 
Pi(E) = —(dE/dx)i on is the rate of the ionization energy losses which, as ever, are assumed 
to be continuous. The symbol (D) denotes a functional describing the "discrete" muon 
energy loss resulting from radiative and photonuclear processes: 

(D(E,x)) = Y / (D(E,x)) k , (2.3) 
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Here da k ' (Ei, E 2 )/dE 2 is the differential cross section for a muon interaction of type k: 



(k = n), and E\ (E2) is the initial (final) muon energy; Nq is the Avogadro number; Z and 
A are the atomic number and atomic weight of the target nucleus. The brackets (. . -)z,a 
indicate an averaging over Z and A. Integrations in Eq. (|2.4j) are performed between the 
limits allowed by the fc-type process kinematics: 

E l,rmn( E ) < E l < E l,m^( E ) > E imm( E ) < E 2 < E% (E) . 



Equation (|2.1| ) does not take into account the muon finite lifetime, which is permissible 
for ultrarelativistic energies and/or for dense enough media^. Moreover, in the equation ( |2.1| ) 
(valid within the "straight-forward" scattering approximation) multiple Coulomb scattering 
and the angular deflection due to inelastic scattering have been ignored. This approximation 
is not so inoffensive but an examination of the problem does not enter the scope of the present 
work. 

A way to include the fluctuation effect due to knock-on electron production by muons will 
be considered later on. An estimation of this effect for DIR has been made by Nishimura [|T2]]. 
According to Ref. |12| the effect leads to an increase of DIR at all depths by approximately 



3% (in the special case of an initial spectrum D (E) oc E 4 ). A reliable analytical method 
for describing the ionization straggling of relativistic muons with incident energies below 



~ 100 GeV in thick absorbers has been suggested recently by Striganov | 29|j , but processes 
others than ionization were not taken into account. 



2 The average decay range of a muon of energy E is given by 

P \ I P 



X d (E) = T^pp/im^c) ~ 6.23 x 10 s g/cm 2 



lg/cm 3 / \lGeV/c 



where m^, r^, and p = ^(E/c) 2 — (m^c) 2 are the muon mass, lifetime, and momentum, respec- 
tively. Clearly Xd(E) is much longer than the muon ionization range \i(E) in a dense medium, 
so the muon decay effect is totally unessential in all instances of interest. 



Let us introduce the macroscopic cross sections by the definition 

/N da z k > A (v J E) \ 
Z k (v,E) = {- J (2.5) 



where 



dv 



Eda z k > A (E,E>) 



dE' 



E'=(l-v)E 

and v is the fraction of energy lost. With help of Eq. ( |2.5|) we rewrite Eq. ( |2.4j) in the more 
convenient form: 

(D{E, x)) k = - v)-^ k (v, E V )D(E V , x) - $ k {v, E)D(E, x)}dv . (2.6) 

Jo 

Here and below 

E) = 6(v^ x (E) - v)6(v - E) , (2.7) 

6{x) is the usual step function, v^ in (E) and u^ iax (£') are the extreme values of v for the 
fc-type process, E v = E / (1 — v), and the function $ k (v, E v ) is defined by Eq. (|2.7|) with the 
substitution E =>• E v . 

At ultrarelativistic energies (E ^> m^c 2 , specifically at E above ~ 10 GeV), we have 
with sufficient accuracy that 

< in (£) = , v k max (E) = 1 , (2.8) 

and hence Q k (v,E) = H k (v,E). Moreover, this may formally be extended to all energies 
considering that radiative and photonuclear losses are inessential to an accuracy of about 



1% at energies under 10 GeV for all media of interest in cosmic ray physics p0|j31[| , and only 
ionization losses are important. Nevertheless, in the following we shall use approximation 
(|2.8| ) only for asymptotic estimations, but we shall assume, if necessary, that v^ in (E) <C 1 
and 1 — v^ aax (E) <C 1 for k = p, b, n in the energy region covered. 



B. Some features of muon-matter interactions at high energies 

A detailed description of the cross sections da^' A (v, E)/dv used in our calculations will 
be presented in a separate publication. For a short review, see Ref. [|3l| . To provide an 
inside into the properties of the radiative processes we have presented in the Appendix 
a very simple parameterization of the -^-dependencies of the (normalized) cross sections 



suggested by van Ginneken p2 |. As one can see from the Appendix, strong energy loss 
fluctuations are more probable in bremsstrahlung. The direct pair production cross section 
goes roughly as 1/v 2 to 1/v 3 over most of the range (v > 0.002). Usually these losses are 
treated as continuous. Nevertheless, as it follows from our estimations, the fluctuation effect 
related to pair production is not exactly negligible and it can prove essential at large depths. 
Hence we will be considering the pair production contribution as discrete, together with 
bremsstrahlung. 

To this must be added that in the limit of complete screening, i.e. for 

(where 7^ is the degree of screening and g m ; n ~ m 2 v /[2E(1 —v)] is the minimum momentum 
transfer), the radiative cross sections are functions of the variable v only (scaling). However 
for values of v which are not too small (namely, at 1 — v <C 1) complete screening occurs only 
at very high energies, E ~ 10 TeV. At lower energies the cross sections grow logarithmically 
with E. 

Unfortunately there is no simple parameterization for the differential cross section of 
the inelastic muon scattering on a nucleus da^ ,A /dv. Moreover, both v- and (especially) 
E'-behavior of the cross section are very model dependent. 

According to the vector-meson-dominance hypothesis da^' A /dv is expressed in terms of 
the total cross section for virtual photon absorption by nucleons and nuclei. A generalized 

j 1 • i i n-\n j\ HA All 1 jii •! j i r j pji 



tions in the diffraction region (low 4- momentum transfers, Q 2 , and large photon energies, v): 
growth with energy of the cross section for nucleon photoabsorption and shadowing effects 
in nuclear photoabsorption. An approximate expression for da^' A / dv has been evaluated in 
the framework of the GVDM by Bezrukov and Bugaev [3q|: 



da n ' (v , E)/dv oc a 1 N{v)F n {v } u)/v . 

Here <j^n(u) is the total cross section for absorption of a real photon of energy v = vE by 
a nucleon. In agreement with accelerator and cosmic-ray experiments p4| . |35 
slowly above v ~ 50 GeV and can be represented approximately as 



<7-yN\y) grows 



Cr 7 7v(z/) 



114.3 + 1.647 In 2 



fib 



,47 GeV, 

The growth of cr 7 7v(V) causes da^' A /dv to depend on the muon energy, E. The function 
F n (v, v) decreases slowly with increasing u, gradually compensating the energy dependence 
of cr 7 7v (a manifestation of the shadowing effect of nucleons inside a target nucleus). Never- 
theless, the logarithmic growth of da^' A /dv quantitatively remains up to E ~ 10 TeV and 
possibly in the asymptotics. The v-dependence of F n (v, v) is rather complicated; for v over 
~ 0.1 it falls off roughly as ln-y with increasing v, thus the fluctuation effect due to this 
process is comparatively large. 

It should be mentioned that the absence of any unitarity constraint allows a very rapid 
(in comparison with the GVDM prediction) increase with energy of the total photoproduc- 
tion cross section as a result of the gluonic structure of the high-energy photon ("minijet 



production mechanism") ||36|| . Although available cosmic-ray data obtained with under- 
ground detectors [Bl] (for v up to ~ 10 TeV) and with EAS arrays p5| (y up to 10 3 — 10 4 



TeV !) do not support this possibility, and what is more, a recent study |37| has shown that 
the calculations of Ref. f36[ strongly overestimate the minijet production contribution at 
v > 10 3 TeV, a significant increase of the photoproduction cross section is still anticipated 



at ultra-high energies. Thus the photonuclear interaction is one of the interesting objects 
for study in future experiments with large underground and underwater detectors. 



III. CONTINUOUS LOSS APPROXIMATION 

Let us at first consider the so-called continuous loss (CL) approximation which is often- 
used for estimations of the CR muon intensity and spectrum under thin enough layers of 
matter (see, e.g., Refs. [0,0] and |38|). It can be obtained from Eq. (|2.1|) by a formal 



expansion of the integrand of expression ( |2.6|) in powers of E v at E v = E, to an accuracy of 



0{v). As a result the functional (|2.3| ) becomes 

(D(E, x)) = £ f\l + E^-)$ k (v, E)D(E, x)vdv . 
, Jo oh, 



Let us define 

L k . (E) 

mm v t 

Clearly b k {E) is the relative partial rate of average energy loss due to the fc-type process, 
and 



b k {E)= / § k {v,E)vdv= / Y, k {v,E)vdv. 

JO Jv k . (E) 



(3(E) = /3i(E) + Ej2h(E) = -{dE/d. 



x 



tot 



is the total rate of energy loss. Thus Eq. ( [2.1|) in the CL approximation takes the form 

—D(E,x) = —[(3(E)D(E,x)}. (3.1) 

Here and below D(E, x) stands for the differential muon spectrum in the CL approximation. 
Similarly I(E, x) and J{x) will stand for integral spectrum and DIR, respectively. 
The solution to Eq. ( ^.1| ) with boundary condition ( |2.2| ) is 

D(E,x) = D (£(E,x) / { i {E ' X)) , (3.2) 



where S(E, x) is the (only) root of the equation X(S,E) = x, and 

r&i dE 



X(E 1 ,E 2 ) = f 

JE- 



Ie 2 (3(E) 

is the average range of a muon with initial energy E\ and final energy E 2 . In other words, 
£(E, x) is the energy which a muon must have at the boundary of the medium in order to 
reach depth x with energy E. It is easily seen that S(E,x) is a monotonically increasing 
function of variables E and x, and the following identities are valid for any x' < x and 
E' > E: 

S(S(E } x'),x- x') = £(E', x - X(E', E)) = S(E } x) . 

It is clear also that S(E,0) = E. 

From Eq. ( p.2[ ) a very nice expression for the integral muon spectrum at depth x can be 
obtained. Let Iq(E) = I(E,0) be the integral spectrum at the boundary, then 

/"oo roc 

I(E,x)= D(E',x)dE' = D (E')dE' = I (S(E,x)) . (3.3) 

JE J£(E,x) 

According to Eq. ( |3.3j ) the expression for DIR, J(x), can be written as 

J{x) = I (£(E t ,x)), (3.4) 

where E t is some detection threshold. It can be argued that the value J(x) is practically 
independent of E t at large depths when E t is sufficiently low (really, when £ t <l TeV). 

In spite of the simplicity and physical transparency of the CL approximation, its range 
of application is fairly restricted. The inadequacy of this approximation is obvious from 
the following simple example. Let the initial spectrum, Dq(E), have a breakoff at some 



energy E mSLX , i.e. D (E) = at E > E max . Then, in accordance with Eqs. (|3.2|) and 
( |3.3| ), D(E,x) = and I(E,x) = at x > X(E mail ,E). It is incorrect, of course, at least 
when X(E milx ,E) < Xd(E). We will demonstrate below, within a simple model, that the CL 



IV. ASYMPTOTIC BEHAVIOR (SPS MODEL) 



We consider here the SPS model mentioned in Sec. |. Let us assume that the functions 
<&k(v,E) and the ionization loss rate, Pi(E), are energy independent, 



$ fc = $fc(u) , A = a = const , 



and the initial spectrum is a power function of energy, 



(4.1) 



D (E) = Dl{E) = CE-W) 



(4.2) 



Moreover muon energies are assumed to be high enough so that conditions Q2.8|) is fulfilled. 

In the SPS model the rate of the average energy loss is simply a + bE, where b = 
bp + K + b n is a constant^, and, therefore, the differential and integral muon spectra in the 
CL approximation are described by 



D SPS (E,x)=D^(E)e 



-■ybx 



l + j£<l-«-) 



1-(7+l) 



(4.3) 



I SPS (E } x) = mE)e 



-■ybx 



1 + iE (l - e 



—bx\ 



(4.4) 



where Iq(E) = 7 _1 Ci?~ 7 is the initial integral spectrum. 

A characteristic property of the SPS model in the CL approximation is a flat (energy 
independent) spectrum (both differential and integral) for E « W e a/i ~ 1 TeV at 
sufficiently large depths (x ^> 1/6), 



3 In reality the quantities (k = p, b, n) and /3j grow with energy logarithmically (or as a power of 



logarithm) up to E ~ 10 TeV (see Refs. |3(],|3l| and p9| ). But, as we have noted in Sec. IIB , it is 
not inconceivable, strictly speaking, that the growth of the relative rate of the average photonuclear 
loss extends at E 3> 10 TeV if a dramatic increase of the total photoproduction cross section with 



and recovery of its original form, 

D SPS (E, x) ~ ^(E)e-^ , 7 SPS (£, x) ~ lS{E)e^ , 

for £7 3> W at all depths. According to Eq. ( |3.4| ), DIR takes the form 

Jsps(x) = U(W{e bx - 1)) , 

independently of the threshold energy E t if E t <C 1^(1 — e~ 6x ). 

Let us consider now the exact solution to Eq. Q2.1|) within the framework of the SPS 
model. Denote 

b 1+n = f 1 §{v)[l - (1 - vy +n ]dv , n = 0, 1, . . . , 
Jo 

and 

£ 7 = 6 7+1 — b ~ = / $(f)(l — v) J vdv , 
Jo 

with $(t;) = + + $n(f)- We shall seek the solution as a series in powers of the 

dimensionless parameter £ = a/{g 1 E): 

D SPS (E,x) = DZ{E)e+* £ (l±ik/ n (x)(-£)» (4-5) 

(here (. . .) n is the Pochhammer symbol). Substituting Eq. (|4.5| ) into Eq. ( |2.1|) , we find that 
the coefficient functions f n (x) satisfy the following recurrence formula: 

fn( x ) + ( b i+n - bry)f n (x) = ng 7 / n _i(i) , /n(0) = <5 n0 • (4.6) 

Integration of Eq. ( |4.6| ) yields 

/„(x) = S n0 + ng 1 I exp[-(6 7+n - b 1 ){x - x')]f n -i(x)dx' . (4.7) 



In particular, for n = and 1 we have from Eq. (|4.7|) fo(x) = 1 and f\(x) = 1 — e~ 0jX . 

By induction, and using the fact that b 7+n — 6 7 < ng^ at n > 1, one can easily verify 
that 

[/l(x)] n < f n (x) < ( Qj x) n 

for all values of x. Therefore the series (|4.5| ) is absolutely and uniformly convergent under 
the condition 

C = = f < 1 , (4.8) 

but it is divergent when £fi(x) > 1. 

It can be shown that the obtained solution reduces to the solution in the CL approxi- 
mation ( |4.$ ) if one sets formally b 1+n = (7 + n)b, for n > 0. A rough fulfillment of these 
equalities at not too large values of n, which is a consequence of a quick growth of the 
electrodynamic cross sections at t) < 1 (see the Appendix), serves as the basis for the ap- 
plicability of the CL approximation. It is obvious, however, that for any n > and 7 > 1 
the exact inequalities b 1+n < (7 + n)b take place, which are satisfied (irrespective of the 
behavior of the function $(t> )), in so far as (1 — v) f > 1 — tv, at any t > 1 and < v < 1. 

Thus the ratio 

r(E,x) = D SPS (E,x)/D SPS (E,x) , 

which is a measure of the fluctuation effect, increases with depth as exp[(7& — 6 7 )x] at 
^ < 1. In other words the CL approximation underestimates the muon intensity at high 
energies. The magnitude of the effect depends critically on the slope of the initial spectrum 
(the remainder 76 — 6 7 quickly increases with 7) and it can be very large. To cite a single 
example, r(E, x) is about 10 at E = 10 TeV and x = 10 km of water equivalent (for standard 
rock), in the case of the vertical spectrum of conventional CR muons from the decay of it 



and K mesons |14]]. It should be noted at the same time that the ratio r(E,x) does not 
necessarily exceed unit at all energies. 

The model under consideration shows that it is impossible to take into account the 
fluctuation effect on muon spectra at large depths as a correction to the CL approximation 
and a reliable method is required. Clearly the exact solution ( |4.5|) by itself is unsuitable for 
calculations at fairly low energies and/or at large depths; it cannot be used, in particular, 
to compute DIR. At the same time the SPS model suggests a starting point for the required 
method: we may, using an ansatz which has the correct asymptotic behavior at high energies, 
construct the solution for the TE applying an appropriate iteration procedure. In the next 
sections we will consider this approach to the problem. 

It will be convenient to specify the asymptotic behavior of the cross sections and ini- 
tial spectrum at high energies as in the SPS model, that is to demand the fulfilment of 
equalities ( |4.1| ) and ( }4.2| ) at energies E ^> E as , where _E as (a conventional bound of the 
asymptotic regime) is a sufficiently large quantity. Then the SPS model will serve as a base 
for asymptotic estimations. It will be recalled that the asymptotic form of the photonuclear 
cross section contribution Q n (v,E) is actually unknown as well as, strictly speaking, the 
high-energy behavior of the initial muon spectrum, Dq(E). Nevertheless, the condition im- 
posed does not restrict generality as long as a concrete value of the bound of the asymptotic 
regime, _E as , is not indicated. Evidently this condition does not play a part in calculation of 
D(E,x) at E < E^ due to the fast decrease with energy of the initial spectrum. 



V. GENERAL CASE: FIRST APPROXIMATION 



Consider the general case. Assuming analyticity of the ratio D(E v ,x)/D (E v ) as a 
function of the variable v at the point v = 0, let us expand this function in a power series 
in v. This yields 



D{E v ,x) = D (E V ) 



1 + X> n 4 



n=l 



[D(E,x)/D {E)\ 



where 



a _"(n-l\E l g 



Then, introducing the definitions 

A n (E) = f 1 $(v, ^(v, E)v n dv , n = 1, 2, . 
Jo 

^(S) = / W, - EM", E v )\dv , 
Jo 



(5.1) 



(5.2) 



with r](v,E) = (1 - v)- 1 D (E V ) / D (E) , we find 



jr A n (E)D (E)d n Du\E) - A(E) 



Ln=l 



£>(£,x) . 



(5.3) 



Due to the fact that the functions £7) depend rather slowly on £7, and the initial 
muon spectrum D (E) is close to a power- low one at high enough energies (vide supra), 
the ratio D(E, x)/D (E) should be asymptotically a relatively slowly varying function of E. 
Thus the derivatives D (E)d n D Q ~ 1 (E)D(E, x) are small. It is obvious also that the integrals 
A n (E) decrease with increasing n. Moreover, due to the specific v -dependence of the cross 



sections (see Sec. |IIB|) , A\(E) ^> A n (E) at n > 1. These simple heuristic considerations 
allow us to use as a first approximation only two leading terms of the expansion ( |5.3|) . In 



this approximation Eq. (2.1) is merely a partial differential one, 



^--8 X (E)^-+K(E) DW(E,x) = 0, 



(5.4) 



where the following notations has been used: 



Pi(E) = (3i(E) + A 1 (E)E 



K{E) = A(E) - [g(E) + l]Ai(E) - fi{E) 



with g(E) + 1 = -EDq 1 {E)D' (E). We will assume subsequently that g(E) is a positive 
definite and nondecreasing function. Clearly g(E) = 7 as £ > E as . 



The function S\(E,x) can be obtained from the equation Xi(£i,E) = x (an analog of 
the energy-range relation), with 



The properties of the function £i(E, x) are completely similar to the ones of above-mentioned 
function £(E,x), but the physical meaning of this quantity is not so obvious. Considering 
that the function j3\{E) is an effective rate of the average energy loss (both continuous and 
discrete) for a given initial muon spectrum, the function £i(E, x) can be interpreted as the 
effective (for the given D (E)) energy which a muon must have at the boundary of the 
medium in order to reach depth x having energy E with a nonzero probability. To refine 
this interpretation let us rewrite Eq. ( |5.5|) in the form which is like the expression for the 
spectrum in the CL approximation ( |3.2| ): 



The solution to Eq. ( |5.4|) can be expressed as 



D {1 \E,x) =D (8i(E,x))exp[-K(E,x)] =V(E,x) , 



(5.5) 



where 




(5.6) 




V(E,x) = D (£ l (E,x)) 



Pi(£i(E,x)) 
HE) 



V(£i(E,x),E), 



(5.7) 



where 



V{E 1 ,E 2 )=exp 



E i q(E')dE' 
e 2 ~MW) 



and 



q(E) = 11(E) + (3[(E) = A(E) - g(E)A 1 (E) + A[(E)E 



(5. 



(5.9) 



Evidently the function q(E) reflects the effect of muon range straggling. It can be demon- 
strated that q(E) > at least for high enough energies. Indeed, substituting Eqs. ( |5.1| ) and 
( |5.2| ) into the right side of Eq. ( |5.9|) yields 

q(E) = [\$(v,E)-[l+g(E)v]7 t (v,E)$(v,E v )}dv 



+ 



o 



\g{E v ) - g{E)]<$>{v, E v ) v (v, E)vdv 



1 „ d$(v,E v ) . ^ , 
E v — \^ E)vdv 



. )r , (5-10) 

io dE v 

The factor [1 + g(E)v]r)(v, E) does not exceed unit[] and decreases fast (tends to zero) 
with increasing v, while the function <&(v,E v ) depends on the second argument, E v , only 
logarithmically. Thus the first integral in Eq. ( |5.10| ) is positive. The second integral is 
nonnegative on the assumption that g(E) is an increasing (or constant) function. The third 
integral is small in comparison with the first one due to the factor rj(v, E)v in the integrand 
and (mainly) to the inequality 

d$(v, E v 



E n , 



0E V 



< $(v,E v ) , 



which takes place even in the absence of the full screening [notice that jz( v ,E v ) < 1 at E 
above ~1 TeV at any v}. Hence the last contribution cannot change the sign of the function 
q(E). 



4 Because the derivative 

^{[1 + g(E)v] V (v, E)} = -g(E v ) {l - ^ + [1 + g^)]^}^, E) 

i • r r\ 1 / r\ 7-l\ -i 



Thus the function q(E) can be interpreted as an effective absorption coefficient dependent 
upon the radiative and photonuclear energy losses, and the function V(£i(E,x), E) should 
be treated as the probability for a muon with energy £\(E, x) at the surface to reach depth 
x with energy E. 

Simple examination shows that A\(E) < b(E). Therefore, £\(E,x) < £(E,x) for all 
values of E and x. Moreover, the remainder £(E,x) — £\(E,x) increases fast with depth 
since 

^[£(E,x)-£ 1 (E,x)] = (3(£(E,x))-p 1 (£ 1 (E,x)) 

= Pi(8) - + K£)£ - Ai^i > , 

and we have taken into account that fli(E) is a nondecreasing function of E after a broad 
minimum at p ~ 300 MeV/c, almost independently of the medium JHJ. It is obvious 
also that the remainder £(E,x) — £\(E,x) increases when the slope of the initial muon 
spectrum grows. The decrease of the minimal muon energy at the surface, necessary in 
order that a muon can reach a given depth with a given energy, is an evident reflection of 
the discreteness of radiative and photonuclear muon energy losses. The function £i(E,x) is 
a useful approximation to estimate this minimal energy within the scope of the approximate 
solution ( |5.5D . 

From Eqs. ( |5.7| - |5.9| ) the following expression for the integral spectrum in the first ap- 



proximation can be obtained: 

I {l) (E,x)= / D {E')V(E',£ x (E',-x))dE' , 

J£i(E,x) 

which is an evident generalization of Eq. ( ^.3[ ) obtained in the CL approximation. It is 
obvious that I^(E,x) < Io(£i(E,x)). In the realistic case, when q(E) is a function slowly 
varying with energy we find 



I^(E,x) ~/ (f 1 (£;,a;))e-^ , 

where g is the average of q(E). 

Consider now the approximate solution ( |5.5|) in the SPS model. It is clear that all 
moments 



A n = f 1 $(«)(1 - vVv n dv = Al 
Jo 



(in particular, Ai = A^ = f? 7 ) and the parameter A = 6 7 are constant in this case. So the 
effective absorption coefficient q = 6 7 — 7£> 7 = g 7 is a positive constant, such that 

^7(7 + 1)A3 < q y < 7 2 A° . 

One can easily show that 



£! PS (E, x) = E[(l + i)e B ^ x - f] and V(£^(E, x), E) = e 



SPS/ 



-q 1 i 



By this means the differential and the integral spectra can be written as 



D$> s (E,x) = D2(E)e-»^ l+£(l-e 



-(7+1) 



(5.11) 



lj$ a (E,x)=l2(E!)e-*r* 1 + ai-e 



(5.12) 



and DIR becomes 



4&G»0 = UiW^* - l))e~** 



for any E t <C W 7 (1 — e e7:r ), where W 7 = o/g 7 . Thus, at large depths, 



^sps^/^sps^) ^ (Qj/bye 



As might be expected, the first two terms of the exact 1/E -expansion ( |4.5| ) coincide with 



,(1) /rn 



r i i 



( |5.11| ) has the correct behavior at high energies at least when £ < 1 (see fl4.8|) ). It should 
be noted also that the approximation ( |5.11| ) is self-consistent. Indeed, one can show that 



d„ 



where 



Z(E,x) 



D l Ss(E,x) 

mm 



Q 1 X\ 



i+e(i-e-^) ' 

Considering that Z(E, x) < 1 at any E and x the series in the right side of Eq ( |5.3|) is always 
uniformly convergent and one may actually cut it off after the 1st term if 

and that is certainly admissible when £fi(x) 1. This is supporting the approximation 
MO. 01 J clS db suitable ansatz. 



VI. GENERAL CASE: ITERATION SCHEME 



Let us now represent the solution to Eq. ( |2.1|) by the following form 



D(E,x) = D {1 \E,x)[l + 5(E,x)} 



(6.1) 



where 5(E,x) is an unknown function ("relative correction"). To derive the equation for 
S(E,x) it is convenient at first to rewrite Eq. (|5.4j ) as 

(6.2) 



— V(E,x) - —\J3i(E)V(E,x)] = [A 1 (E)u(E,x)-A(E)]V(E,x) 



where 



w(E, x) 



ME) 



(6.3) 



h(E) = n(E) + M^±M = A{E) + m+mm - m) . (6.4) 

E E 



In order to derive Eq. ( |6.2| ) we have used Eqs. (|5.5| ) and ( |5.6|) . Direct substitution of Eq. (|6] 
into Eq. ( |2.1| ) , in view of Eq. ( |6.2| ) , then gives 

Ci6(E, x) = f 1 $(v, E v ) {Q(E, x; v) [1 + 5(E V , x)} 
Jo 

- [l + u(E,x)v)[l + S(E,x))}r]{v,E)dv , (6.5) 



where the differential operator 



was introduced, and we have defined 

Q(E, x- v) = ^j§L D "f F { f;^ ex P [/C(E, x) - K(E Vi x)] . (6.6) 
D (E V ) D Q (£ x (E,x)) 

Clearly 5(E, 0) = 0. We shall seek the solution to Eq. ( |6 . 5| ) using a procedure of successive 
approximations. 

Let us note initially that the function S(E, x) follows a C2(x) / E 2 -dependence as E ^ E^, 
where C2(x) is independent of energy. It is a straight corollary of the coincidence of the first 
two terms in the 1/E -expansions for the approximate solution ( |5.5|) and the exact SPS 
solution ( |4.5|) . Therefore 

B(E, X] v) = 5(E V , x)-(l- v) 2 5(E, x) oc (1 - v) 2 v/E 3 

as E ^> E^, i.e. the function Q(E,x;v) is small in absolute value by comparison with 
S(E, x). We assume that at all energies the term with the factor Q(E, x; v) can be neglected 
in the integrand of the right side of Eq. ( |6.5| ) as a first approximation. Thus the equation 
for the correction function in second approximation becomes 

r r> n f T7i \T rf2"l / 7-1 \ n f tti \ cCZ) / ttt n\ r\ I ^ f-r\ 



where 



Ri(E, x) = $0, E v ) {fl(E, x\ v)(l - v) 1 - [1 + u(E, x)v]}rj(v, E)dv . (6. 



Solving Eq. ( |6.7| ) yields 



<5 (2) (£,x) / 0X1) 
10 



R 2 (£i(E, x — x"),x")dx" 



Ro(£i(E, x — x'),x')dx' 



Ei{E,x) 



exp 



W*) R^E^x-HE.E")) ^, 



R (E>x-X^ EjldE ^ 



where Si(E, x) is the only root of the equation Aj(£j, E) = x, and 

\{E\,E 2 ) 



Ie 2 /3i(E) 

is the ionization range of a muon with initial energy E\ and final energy Ei (hence £i(E, x) — 
E is simply the energy lost due to ionization). 

Let us consider one evident consequence of Eq. (|6.9|) . Clearly R 2 (E,x) < Ro(E,x) for 
all values of the arguments. Substituting this inequality into Eq. ( |6.9|) , and integrating over 
E then gives 



with 



exp[K 2 (E,x)} < l + 5 {2 \E,x) < exp[K {E,x)} 



Ki(E,x)= f X Ri(£i(E,x - x'),x')dx' 
Jo 



(6.10) 



The exponential factors in (|6.10| ) can be treated as the lower and upper limits for the 
correction to the "survival probability" V(£\(E, x), E) so long as 

riqi^E^- x')) - R (£ t (E,x- x'),x')]dx' > . 
Jo 

In order to build an equation for calculation of the correction function in the Z-th ap- 



Cs(x)/E 3 with an .^-independent function c^(x), as it can be easily verified using Eqs. (|6.5| ) 
and ( |6.7| ). Therefore, in the next approximation we may put approximately 

S(E v ,x) -5 {2) {E v ,x) ~ (1 -v) 3 [5{E,x) -5 (2) {E,x)} . 

Repeating the consideration we find by induction that 5(E,x) — 5®(E,x) — > ci(x)/E l as 
E 3> -E as , and thus we put 



6(E v ,x) -6 {l) (E v ,x) ~ (1 -v) l+1 [6(E,x) -5 {l) (E,x)} . 



(6.11) 



Let us define 



Qi(E, x) = 5 {l \E, x) - 5 {l - 1] (E, x) , I > 2 , 



(6.12) 



with 5^(E, x) = by definition. From Eq. ( |6.5[ ), using Eq. (|6.12| ) and Eq. ( |6.11| ) we obtain 
the recursion chain of equations for the functions &i(E, x) : 



- R t {E, x)]Gi(E, x) = $ti-i{E, x), l>3 



(6.13) 



where 



?fl l (E,x)= / ^(v,E v )Q(E,x',v)[Qi(E v ,x) - (l-v) l ei(E,x)]r](v,E)dv 
Jo 

The solution to Eq. (|6.13|) is given by 

px r px 

Qi(E,x)= / exp / RAEAE, x — x"), x")dx 

Jo \_Jx> 



(6.14) 



?Ri-i(£i(E, x — x ), x')dx' 



£i(E,x) 



cxp 



Ri(E", x - Xi(E, E^JI^h 



jw) dE - ( ° 5) 



To verify the convergency of the iteration procedure consider firstly the behavior of the 
function Ri(E,x) at I 3> 1. Due to the factor (1 — v) 1 in the first term of the integrand of 
Eq. (Bit) and the properties of the macroscopic cross sections (see Sec. 11 Bl), onlv the region 



of small values of v is important in this case. So at I ^> 1 the function Q(E, x; v) can be 
estimated as 



n(E,x;v) ~ tt(E,x;0) + v 
or, considering the definitions ( |6.3|) , (|6.4f ), and (|6.6|) 



dQ(E, X] v) 
dv 



J D = 



x; u) ~ 1 + ac)u ~ 1 . 



(6.16) 



Thus 



R t (E, x) ^ - / $(v, £„) [1 + [1 - (1 - 

•/ o 

at / ^> 1 and, therefore, Ri(E,x) < and \Ri(E,x)\ increases indefinitely with I. Clearly 
the exponential factor 



exp 



f RA£AE,x-x"),x")dx" 

Jx' 



in the integrand of Eq. (|6.15|) diminishes fast with increasing I. 

On the other hand, in view of the fact that Qi(E,x) oc E~ l at energies high enough, we 
can write Qi(E v , x) — (1 — v) l Qi(E, x) = v(l — v) l Fi(E, x; v), where Fi(E, x; v) is a function 
which can be estimated at v <C 1 by 

,ae,(£,x) 



F,(S,z;0) = -E- 



<9£ 



lQi{E,x) oc . 



Hence, using Eq. ( |6.16| ), the integral (|6.14j) can be estimated at I ^> 1 as 

Mt(E,x) ~ Fj(£,x;0) [ $(v,E v )[l + v(E,x)v]r)(v, E)(l - v) l vdv . 

Jo 

Thus x) and ^i{£i{E,x — x'),x') are positive and decrease when / increases, if the 

function Fi(E, x; 0) is bounded in magnitude as I — > oo or even if x; 0)| increases with 

/ not too fast. This can be verified by induction. 



The foregoing proves that Qi(E, x) —>■ as I — > oo for any depths at least for high enough 
energies. Due to the correct asymptotic behavior of the functions Qi(E, x) at all values of / 
this indicates that the iteration procedure converges, that is 

6(E,x) = \im5 (l) (E,x) . 

I— >oo 

A more cumbersome analysis and numerical verifications demonstrate that this statements 
is true at all energies under quite general assumptions on the energy dependence of the 
rate of the continuous energy loss, the macroscopic cross sections, and the initial spectrum; 
specifically if the functions (3i(E) and E) increase monotonically and sufficiently slowly, 
while Dq(E) decreases with energy so that g(E) is a slightly varying function of energy. As 
it follows from numerical estimations, the convergency rate is very good, and usually only a 
few iterations are needed to reach an accuracy of the order of 1% at depths up to ~ 20 km 
of water equivalent for muon energies above ~ 1 GeV. 

By way of illustration let us again direct our attention to the SPS model and consider the 
second approximation correction function 5gp S (E,x). Under condition (|4.8|) we can write 
the exact expression for the correction function, using definition ( |6.1| ) and Eqs. ( |4.5| ) and 



^■^I ^FalP '- '' (6 - 17) 

Therefore at £ < 1 the leading asymptotic term has the form 



c 7 e 



a + a x e-^ x + a^ 2 ^ - (a Q + a 1 + a 2 )e- {2 - e)e ~< x f , (6.18 



-2 



where c 7 = (7 + 1)(7 + 2)/2, a = 1/(2 - e), a 1 = -2/(1 - e), a 2 = -1/e, and e = AJ/A?. 
The expression for Ri(E,x) in the SPS is given by 

R t (E,x) = W l (Z{E,x)) (6.19) 



(i- v y 



R]{z) = J o $(„) | ^_ zv > i+1 - [1 + (7 + l)zv\\ (1 - vfdv . 
In particular, at 2 < 1 we have 

Rl{z) ~ c 7 A]z 2 , Rl{z) ~ -(2 - e)A7 . (6.20) 



One can easily show, using Eqs. ( 6.19Q and (|6.20|) that the asymptotic behavior of the exact 



correction function ( 6.18| ) is reproduced, as was to be expected, by the correction function in 



the second approximation (see Eq. ( p.9|) ). It is important that the approximate correction 

(2\ 

8$p S (E, x) is definite and bounded at all values of E and x, in contrast with the exact 
expression ( p.l7| ) given by a series which converges only under condition ( |4.8| ). This can be 
seen from the constraint ( |6.10 ). Indeed, a little manipulation, and taking into account that 



yields 



dZ(E, x)/dE < and dZ(E,x)/dx > 



00 fry _|_ 1 1 rx 

K (E } x) = Y W , ,n AZ / Z n (E + ax', x - x')dx' 
^2 nl J o 



< xR 1 (Z(E,x))Z(E + ax,x)/Z(E,x) 

< xR2 ( ] 1 + ^ - = K™ X (E, x) . 



Hence 



< 5§ s (E,x) < exp[Kr x (E,x)} - 1 . 

and, therefore, the survival probability calculated in the second approximation, does not 
exceed the factor exp{ — [g 7 x — K™ a *(E,x)}}. It is obvious that the function K™ ax (E,x) 

(2) 

and, therefore, the correction 5gp S (E,x) are both bounded at all finite values of E and 

(2) 

x. A significant consequence resides in the fact that <5 SPS <C 1 at £ <C 1 for any x since 
^-max ^ c ^ e £^/(i <^ i i n other words, the first approximation solution ( |5.11| ) is correct 



Analogous statements can be also proved in the general case: the first approximation 
solution ( |5.5[ ) is practically exact at all depths when E 3> @i(E)/Ai(E). 



VII. SUMMARY AND OUTLOOK 



The method described enables us to calculate with a controlled precision the differential 
energy spectra of CR muons after propagation through thick layers of matter. It is appro- 
priate for any depths when muon energies are high enough and provides a way for including 
the real (non-power-law) initial muon spectrum and the energy variation of the continuous 
and discrete muon energy losses, with only the formal (but rather natural) requirement for 
the asymptotic behavior of the initial spectrum and the cross sections. A computer imple- 
mentation of the method is fully straightforward and the required CPU time is small. This 
enables to use the method in on-line processing for underground/underwater experiments. 
It is important that the useful notion of the minimal muon energy at the boundary £(E, x), 
which has a physical sense when the range straggling is negligible (i.e. for small depths), 
has an analog (£i(E, x)) in the case of arbitrary large fluctuations (i.e. for arbitrary depths) 
when the first approximation solution V(E,x) is available (high muon energies, specifically 
E > a few TeV). 

We intend in the near future to give a detailed numerical illustration of the convergency 
of the iteration procedure besides calculating results on the CR muon energy spectra (dif- 
ferential and integral) and DIR at large depths for different types of rocks and water with 
varying charm production models etc. In recent years a rather representative array of data 
on DIR in rock and (to a lesser extent) in water has been accumulated by many experiments. 
One should systematize all these data and compare them against theoretical predictions. 

We also look forward to analyze future underground and underwater experiments to 
possibly throw light on the prompt muon problem and to derive information about super- 
high-energy muon interactions with nuclei, primarily about the rather poorly studied pho- 
tonuclear interaction. 
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APPENDIX: 

van Ginneken's parameterization of the v-distributions 
for radiative processes 

We give here a parameterization of the normalized cross sections fk(v, E) = a^ l dak/ dv 
(where dak/dv = da^' A /dv) for pair production and bremsstrahlung as a function of the 
fractional energy loss v suggested in Ref. ||32| . The formulas presented below are valid at 
muon energies from ~ 100 GeV up to 30 TeV and enable one to estimate the comparative 
probabilities of "soft" (v <C 1) and "hard" (v ~ "W max ~ 1) losses in the radiative processes. 

• Direct pair production off both nuclear and electron targets • 

f( n ' e \v,E) oc const, 5m e /E< v <25m e /E, 

oc v~ l , 25m e /E < v < v\ (if E > 25m e /v 1 ) , 

oc v ~ 2 , V\ < v < i>2 (if E > 25m e /t> 2 ) , 

oc v ~ 3 , l>2 < v < 1 , 

where v± = 0.002 and v 2 = 0.02 . For very small v (up to the kinematic limit t> ^ in (E) = 
i.m e /E) da p (v, E)/dv oc a (vE) ln(l/t>)i7 _1 , where (Tq{v) is the total cross section for pair 
production by a photon of energy v (Kel'ner's approximation). Thus da p /dv follows roughly 



Below v ps 5m e , ao(v>) remains small (< 0.05<To(oo)), and it increases roughly linearly until 
v ~ 25m e , where <7q(v) ~ 0.5<Jo(oo). 

• Bremsstrahlung contribution off a nuclear target • 

ft\v,E) oc u- 1 , < in (^)< v <0.03, 

oc u c »( B ) , 0.03 < v < v£(E) , 

oc (1 - v) C "W ? v «(£^ < v < 1 , 

where < in (£) = 0.001/E, = (1 + 4.5/v 7 ^)" 1 , C n {E) = 1.39 - 0.024 In E, and 

C' n {E) = 1.32 — 0.12 \nE . Only for values of v above 0.995 the parameterization is not very 
reliable, but it is not important in practice. 

• Bremsstrahlung contribution off atomic electrons • 

ft\v,E) oc v~™\ <;„(#)< v <0.05, 

OC V e {E) , 0.05 < v < v e b {E) , 

oc (v b max (E) - ^) 1/2 , vt(E) < v < < ax (£) , 

where < ax (£) = (1 + 10.92/.B)- 1 , v e b (E) = <(E)< M (£) , and C e (E) = 1.50 - 0.03 In E . 

The proportionality factors which was omitted in the above parameterization can be 
obtained by continuity and normalization. A slight Z-dependence in the ^-distributions has 
been ignored. The energy E and the electron mass m e have been expressed in GeV. 



REFERENCES 



[1] V. S. Berezinsky and G. T. Zatsepin, Usp. Fiz. Nauk 122, 3 (1977) [Sov. Phys. Usp. 
20, 361 (1977)]. 

[2] V. S. Berezinsky, S. V. Bulanov, V. A. Dogiel, V. L. Ginzburg, and V. S. Ptuskin, 
Astrophysics of Cosmic Rays, edited by V. L. Ginzburg (North- Holland, Elsevier Science 
Publishers B. V., 1990). 

[3] A. Roberts, Rev. Mod. Phys. 64, 259 (1992). 

[4] BAIKAL Collaboration, S. D. Alatin et at, in TAUP'91, Proceedings of the 2nd In- 
ternational Workshop on Theoretical and Phenomenological Aspects of Underground 
Physics, Toledo, Spain, 1991, edited by A. Morales, J. Morales, and J. A. Villar [Nucl. 
Phys. B (Proc. Suppl.) 28A, 491 (1992)]. 

[5] L. K. Resvanis (for the NESTOR Collaboration), in Proceedings of the Workshop on 
High Energy Neutrino Astrophysics, Honolulu, Hawaii, 1992, edited by V. J. Stenger et 
al. (World Scientific, 1992), p. 325. 

[6] AMANDA Collaboration, S. W. Barwick et al., in Proceedings of the 22nd International 
Cosmic Ray Conference, Dublin, Ireland, 1991, edited by M. Cawley et al. (The Dublin 
Institute for Advanced Studies, Dublin, 1991), Vol. 5, p. 658. 

[7] H. Sobel, in TAUP'91 §, p. 365. 

[8] G. T. Zatsepin and E. D. Mikhalchi, J. Phys. Soc. Japan, 17, Suppl. A - III, 356 (1962); 
E. D. Mikhalchi and G. T. Zatsepin, in Proceedings of the 9th International Cosmic Ray 
Conference, London, UK, 1965, (The Institute of Physics and The Physical Society, 
London, 1965), Vol. 2, p. 997. 

[91 V. I. Gurentsov. G. T. ZatseDin. and E. D. Mikhalchi. Yad. Fiz. 23. 1001 fl976) fSov. J. 



Nucl. Phys. 23, 527 (1976)]; V. I. Gurentsov, Institute for Nuclear Research (Academy 
of Sciences of the USSR) Report No. II - 0379, Moscow, 1984 (unpublished). 

[10] J. Nishimura, in Proceedings of the 8th International Cosmic Ray Conference, Jaipur, 
India, 1963, edited by R. Daniels et al. (Commercial Printing Press, Ltd., Bombay, 
India, 1964-5), Vol. 6, p. 224; in Handbuch der Physik (Springer - Verlag, Heidelberg, 
1965), Bd., 46/2. 

[11] S. Hayakawa, J. Nishimura, and Y. Yamamoto, Prog. Theor. Phys. Suppl. 32, 104 
(1964). 

[12] J. Nishimura, in Proceedings of the 9th International Cosmic Ray Conference ||, Vol. 
2, p. 1003. 

[13] A. Misaki and J. Nishimura, Uchuusen-Kenkyuu, 21, 250 (1976). 

[14] E. V. Bugaev, V. A. Naumov, and S. I. Sinegovsky, Institute for Nuclear Research 
(Academy of Sciences of the USSR) Report No. II - 0347, Moscow, 1984 (unpublished); 
Yad. Fiz. 41, 383 (1985) [Sov. J. Nucl. Phys. 41, 245 (1985)]. 

[15] E. V. Bugaev, V. A. Naumov, and S. I. Sinegovsky, Izv. Akad. Nauk SSSR. Ser. Fiz. 
49, 1389 (1985) [Proc. of the Acad, of Sci. of the USSR. Phys. Ser. 49, 146 (1985)]; see 
also Ref. . 



[16] K. Kobayakawa, Nuovo Cimento B 47, 156 (1967); E. Kiraly, P. Kiraly, and J. L. 
Osborn, J. Phys. A 5, 444 (1972); Y. Minorikawa, T. Kitamura, and K. Kobayakawa, 
Nuovo Cimento C 4, 471 (1981); D. P. Bhattacharyya, Nuovo Cimento C 9, 404 (1986); 
O. C. Allkofer and D. P. Bhattacharyya, Phys. Rev. D 34, 1368 (1986); V. N. Bakatanov, 
Yu. F. Novosel'tsev, R. V. Novosel'tseva, A. M. Semenov, and A. E. Chudakov, Yad. 
Fiz. 55, 2107 (1992) [Sov. J. Nucl. Phys. 55, 1169 (1992)]. 



[17] S. Miyake, V. S. Narasimham, and P. V. Ramana Murthy, Nuovo Cimento 32, 1524 
(1964); H. Oda and T. Murayama, J. Phys. Soc. Japan 20, 1549 (1965); ERPM Group, 
B. S. Meyer et al, Phys. Rev. D 1, 2229 (1970). O. C. Allkofer, D. P. Bhattacharyya, 
and P. Pal, Nuovo Cimento C 9, 421 (1986); 

[18] A. A. Lagutin, A. G. Prokopets, A. K. Konopelko, and S. I. Sinegovsky, in Proceedings 
of 22nd International Cosmic Ray Conference 0, Vol. 2, p. 752. 

[19] P. J. Hayman, N. S. Palmer, and A. W. Wolfendale, Proc. Roy. Soc. (London) A275, 
391 (1963); J. L. Osborn, A. W. Wolfendale, and E. C. M. Young, J. Phys. A 1, 55 
(1968); L. Bergamasco and P. Picchi, Nuovo Cimento B 3, 134 (1971); Yu. N. Vavilov, 
Yu. A. Trubkin, and V. M. Fedorov, Yad. Fiz. 18, 884 (1974) [Sov. J. Nucl. Phys. 18, 
434 (1974)]; J. N. Capdevielle, D. P. Bhattacharyya, B. S. Acharya, J. Procureur, M. 
F. Bourdeau, and P. Pal, J. Phys. G 11, 565 (1985); V. A. Kudryavtsev, Institute for 
Nuclear Research (Academy of Sciences of the USSR) Report No. II - 0529, Moscow, 
1987 (unpublished); H. Bilokon, A. Castellina, B. d'Ettorre Piazzoli, G. Mannocchi, P. 
Picchi, and S. Vernetto, Nucl. Instrum. Methods A303, 381 (1991); K. Mitsui, Phys. 
Rev. D 45, 3051 (1992). 

[20] V. I. Gurentsov, Institute for Nuclear Research (Academy of Sciences of the USSR) 
Report No. IT - 0380, Moscow, 1984 (unpublished). 

[21] N. Takahashi, H. Kujirai, A. Adachi, N. Ogita, and A. Misaki, Uchuusen-Kenkyuu, 28, 
120 (1984). 

[22] M. G. K. Menon and P. V. Ramana Murthy, in Progress in Elementary Particle and 
Cosmic Ray Physics, edited by J. G. Wilson and S. A. Wouthuysen (North-Holland 
Publishing Company, Amsterdam, 1967), Vol. 9, p. 163. 

[23] E. V. Bugaev, Yu. D. Kotov, and I. L. Rozental', Cosmic Muons and Neutrinos (Atom- 



izdat, Moscow, USSR, 1970) (In Russian). 

[24] Yu. N. Andreyev, A. E. Chudakov, V. I. Gurentsov, and I. M. Kogai, in Proceedings 
of the 21st International Cosmic Ray Conference, Adelaide, Australia, 1990, edited by 
R. J. Protheroe (Department of Physics and Mathematical Physics, The University of 
Adelaide, Graphic Services, Northfield, South Australia, 1990), Vol. 4, p. 301. 

[25] H. Inazawa, K. Kobayakawa, and T. Kitamura, Nuovo Cimento C 9, 382 (1986); K. 
Mitsui, Y. Minorikawa, and H. Komori, ibid 9, 995 (1986); L. V. Volkova, W. Fulgione, 
P. Galeotti, and O. Saavedra, ibid 10, 465 (1987). 

[26] E. V. Bugaev, V. A. Naumov, S. I. Sinegovsky, and E. S. Zaslavskaya, Nuovo Cimento 
C 12, 41 (1989); Izv. Akad. Nauk SSSR. Ser. Fiz. 53, 342 (1989) [Proc. of the Acad, of 
Sci. of the USSR. Phys. Ser. 53, 135 (1989)], and references therein. 

[27] R. P. Kokoulin and A. A. Petrukhin, in Proceedings of the 22nd International Cosmic 
Ray Conference ||, Vol. 4, p. 536. 

[28] I. L. RozentaF and V. N. Strel'tsov Zh. Eksp. Teor. Fiz. 35, 1440 (1958) [Sov. Phys. 
JETP 35, 1007 (1959)]. 

[29] S. I. Striganov, Nucl. Instrum. Methods A322, 225 (1992). 

[30] Particle Data Group, K. Hikasa et al, Phys. Rev. D 45, SI, (1992). 

[31] W. Lohmann, R. Kopp, and R. Voss, CERN Yellow Report 85 - 03, Geneva, 1985 
(unpublished). 

[32] A. van Ginneken, Nucl. Instrum. Methods A251, 21 (1986). 

[33] L. B. Bezrukov and E. V. Bugaev, Yad. Fiz. 32, 1636 (1980) [Sov. J. Nucl. Phys. 32, 
847 (1980)]; ibid. 33, 1195 (1981) [33, 635 (1981)]. 



[34] See G. T. Zatsepin, E. V. Korol'kova, V. A. Kudryavtsev, A. S. Mal'gin, and O. G. 
Ryazhskaya, Yad. Fiz. 49, 426 (1989) [Sov. J. Nucl. Phys. 49, 266 (1989)]; V. N. 
Bakatanov, A. E. Chudakov, Yu. F. Novosel'tsev, R. V. Novosel'tseva, A. M. Semenov, 
and Yu. V. Sten'kin, in Proceedings of the 21st International Cosmic Ray Conference 



24j| , Vol. 9, p. 375; and references therein. 



[35] See D. Dumora, J. Procureur, and J. N. Stamenov, J. Phys. G 18, 1839 (1992), and 
references therein. 

[36] M. Dress and F. Halzen, Phys. Rev. Lett. 61, 275 (1988); T. K. Gaisser, F. Halzen, T. 
Stanev, and E. Zas, Phys. Lett. B 243, 444 (1990); R. Gandhi, I. Sarcevic, A. Burrows, 
L. Durand, and H. Pi, Phys. Rev. D 42, 263 (1990); R. Gandhi and I. Sarcevic, ibid 44, 
10 (1991); M. Dress and R. M. Godbole, Phys. Rev. Lett. 67, 1189 (1991). 

[37] J. C. Collins and G. A. Ladinsky, Phys. Rev. D 43, 2847 (1992); R. S. Fletcher, T. K. 
Gaisser, and F. Halzen, ibid 45, 377 (1992), and 45, 3279 [Errata] (1992); J. R. Forshaw 
and P. N. Harriman, ibid 46, 3778 (1992). 

[38] J. N. Crookes and B. C. Rastin, Nucl. Phys. B58, 93 (1973); A. I. Barbouti and B. C. 
Rastin, J. Phys. G 9, 1577 (1983); I. W. Rogers and M. Tristam, ibid 10, 983 (1984). 

[39] L. B. Bezrukov and E. V. Bugaev, in Proceedings of the 17th International Cosmic 
Ray Conference, Paris, France, July 1981 (Section d'Astrophysique, Centre d'Etudes 
Nucleaires de Saclay, Gif-sur-Yvette, Cedex, 1981), Vol. 7, p. 102. 



